
--****
--  author  : jiri babor
--  email   : jbabor@paradise.net.nz
--  project : basic matrix operations
--  date    : 99-09-06
--  version : 1.00

--  (1) Yes, I am aware of Archive routines by Arthur Adamson and
--      J.F.Soucaille, but I prefer to roll my own...
--  (2) I am very much interested in better and/or faster routines,
--      but I am also quite immune to the type check diarrhea a lot
--      of Euphorians seem to be suffering from, so don't even try
--      to inflict it on me. Thanks.
namespace mat

include std/math.e



constant
	bound=1.0e-12,
	true=1,
    false=0
    
public integer error

--~ function abs(atom x)
	--~ if x<0 then return -x end if
    --~ return x
--~ end function

--~ function odd(integer i
	--~ )
	--~ -- return true (1) if i is odd, else false (0)

	--~ return and_bits(i,1)
--~ end function

--**
-- return filled matrix of chosen dimensions
-- initializing fill value default to 0
public function mat(
	integer number_of_rows,
	integer number_of_columns,
	object fill=0	
	)

    return repeat(repeat(fill,number_of_columns),number_of_rows)
end function

--**
-- return identity matrix of given order
public function identity(integer order	-- dimension of square matrix
	)
	
	sequence s

	s=mat(order,order,0)
	for i=1 to order do
    	s[i][i]=1
	end for
    return s
end function

--**
-- return transpose of m
public function transpose(sequence m
	)
		sequence s
	integer r,c

	r=length(m)
	c=length(m[1])
	s=mat(c,r,0)
	for i=1 to c do
    	for j=1 to r do
    		s[i][j]=m[j][i]
		end for
	end for
    return s
end function

--**
-- return a*b, matrix product
public function multiply(sequence a, sequence b)
		sequence s,t
	integer m,n,p

	m=length(a)
	n=length(b)
	p=length(b[1])
	s=mat(m,p,0)
	for i=1 to m do
    	for j=1 to p do
    		t=s[i]
    		for k=1 to n do
    			t[j]+=a[i][k]*b[k][j]
			end for
			s[i]=t
		end for
	end for
    return s
end function

--**
-- return trace of square matrix
-- mtrace name chosen to avoid collision with built-in trace!
public function mtrace(sequence s
    )
			atom t
	t=s[1][1]
	for i=2 to length(s) do
    	t+=s[i][i]
	end for
    return t
end function

--**
-- return submatrix as result of elimination
public function submatrix(
	sequence s,		-- original matrix
	integer i,		-- row to be eliminated
	integer j		-- column to be eliminated
	)
	
	integer m,n

    m=length(s)		-- number of rows
    n=length(s[1])	-- number of columns
    s=s[1..i-1] & s[i+1..m]
    for row=1 to m-1 do
    	s[row]=s[row][1..j-1] & s[row][j+1..n]
	end for
    return s
end function

--**
-- return determinant of square matrix
public function det(sequence m
	)
	
	atom d
	integer n,s

	n=length(m)
	if n=1 then return m[1][1] end if
	s=1
	d=0
	for i=1 to n do
		d+=s*m[1][i]*det(submatrix(m,1,i))
		s=-s
	end for
    return d
end function

--**
-- return inverse of sequence a
-- if global error flag is set then returned value is meaningless
public function inverse(sequence a	-- square nonsingular matrix
	)
		
	sequence b
	integer i,size
	atom factor,temp

    size=length(a)
    b=identity(size)
    for j=1 to size do
    	i=j
        while abs(a[i][j]) < bound do
        	if i=size then
                error=true
        		return b
    		end if
    		i+=1
    	end while
		for k=1 to size do
    		temp=a[i][k]
    		a[i][k]=a[j][k]
    		a[j][k]=temp
		    temp=b[i][k]
    		b[i][k]=b[j][k]
    		b[j][k]=temp
		end for
		factor=1/a[j][j]
		for k=1 to size do
    		a[j][k] *= factor
    		b[j][k] *= factor
		end for
		for m=1 to size do
    		if m != j then
				factor=-a[m][j]
				for k=1 to size do
    				a[m][k] += factor*a[j][k]
    				b[m][k] += factor*b[j][k]
				end for
			end if
		end for
	end for
    error=false
    return b
    end function
    
-- attempt at formatted matrix print
-- output file, 1 prints to screen
-- s matrix to be shown
--mode the printf mode string checks for integer or atom or supply your own
public procedure mpf(      
    integer f,      
    sequence s,      
    object mode=0   )
    atom x

    for r=1 to length(s) do
        for c=1 to length(s[1]) do
            x=s[r][c]
            if mode then
		printf(f,mode,x)
            elsif integer(x) then
                printf(f,"%12d",x)
            elsif abs(x)<1.0e6 and abs(x)>1 then
                printf(f,"%12.3f",x)
            else
                printf(f,"%12.3e",x)
            end if
        end for
        puts(f,'\n')
    end for
    puts(f,'\n')
end procedure

--ne1 chg atom to object for fill type, add tests, 
--rearrange docs add mode to mpf

--add some transpose and ravel things
--add more dimensions z,t, n
--define some enum to access x & y or i & j then z or k
--allow row or col major access to matchup w/c arrays or other lang

